home *** CD-ROM | disk | FTP | other *** search
/ Libris Britannia 4 / science library(b).zip / science library(b) / SCIENTIF / H381.ZIP / GSRC208A.ZIP / DDOT.C < prev    next >
C/C++ Source or Header  |  1992-05-01  |  2KB  |  85 lines
  1.  
  2. double ddot( n, dx, incx, dy, incy )
  3.  
  4. double *dx, *dy;
  5. int n, incx, incy;
  6.  
  7. /*
  8.    Purpose : Inner product dx . dy
  9.  
  10.  
  11.    --- Input ---
  12.  
  13.    n    : number of elements in input vector(s)
  14.    dx   : double vector with n+1 elements, dx[0] is not used
  15.    incx : storage spacing between elements of dx
  16.    dy   : double vector with n+1 elements, dy[0] is not used
  17.    incy : storage spacing between elements of dy
  18.  
  19.  
  20.    --- Output ---
  21.  
  22.    ddot : dot product dx . dy, 0 if n <= 0
  23.  
  24.  
  25.    ddot = sum for i = 0 to n-1 of
  26.    dx[lx+i*incx] * dy[ly+i*incy] where lx = 1 if
  27.    incx >= 0, else lx = (-incx)*(n-1)+1, and ly
  28.    is defined in a similar way using incy.
  29.  
  30. */
  31.  
  32. {
  33.    double dotprod;
  34.    int ix, iy, i, m;
  35.  
  36.    dotprod = 0.;
  37.    if ( n <= 0 )
  38.       return dotprod;
  39.  
  40. /* Code for unequal or nonpositive increments.  */
  41.  
  42.    if ( incx != incy || incx < 1 ) {
  43.       ix = 1;
  44.       iy = 1;
  45.       if ( incx < 0 )
  46.          ix = ( -n + 1 ) * incx + 1;
  47.       if ( incy < 0 )
  48.          iy = ( -n + 1 ) * incy + 1;
  49.       for ( i = 1 ; i <= n ; i++ ) {
  50.          dotprod = dotprod + dx[ix] * dy[iy];
  51.          ix = ix + incx;
  52.          iy = iy + incy;
  53.       }
  54.       return dotprod;
  55.    }
  56.  
  57. /* Code for both increments equal to 1.  */
  58.  
  59. /* Clean-up loop so remaining vector length is a multiple of 5.  */
  60.  
  61.    if ( incx == 1 ) {
  62.       m = n % 5;
  63.       if ( m != 0 ) {
  64.          for ( i = 1 ; i <= m ; i++ )
  65.             dotprod = dotprod + dx[i] * dy[i];
  66.          if ( n < 5 )
  67.             return dotprod;
  68.       }
  69.       for ( i = m + 1 ; i <= n ; i = i + 5 )
  70.          dotprod = dotprod + dx[i] * dy[i] + dx[i+1] * dy[i+1] +
  71.                       dx[i+2] * dy[i+2] + dx[i+3] * dy[i+3] +
  72.                       dx[i+4] * dy[i+4];
  73.       return dotprod;
  74.    }
  75.  
  76. /* Code for positive equal nonunit increments.   */
  77.  
  78.    for ( i = 1 ; i <= n * incx ; i = i + incx )
  79.       dotprod = dotprod + dx[i] * dy[i];
  80.    return dotprod;
  81.  
  82. }
  83.  
  84.  
  85.